
#include "Transform.h"
#include "Geometry.h"
#include "Sampler.h"
#include "PNTriangle.h"
#include "Image.h"
#include "ImageIO.h"
#include "Mesh.h"
#include "MeshIO.h"

#include <iostream>


using namespace gk;
using namespace std;

static gk::Sampler sampler;
static int colorB;
static int colorG;
static int colorR;
static int nbDehors;
static int nbPixelEcrit;
static int nbPixelNonEcrit;
static float taillePixel = 0.015;
static int nbTriangle = 0;

/// Eclairages
    vector<Point> sourcesLumiere;

    Vector intensiteAmbiante(200,0,100);
    Vector intensiteDiffuse(255,255,255);
    Vector intensiteSpeculaire(255,255,255);
    double coeffSpeculaire = 2;

    float fMatiere = 0.9 / M_PI;

    int largeur = 1024;
    int hauteur = 1024;
    double maxi = 0;

///Prototypes

// Optimisation
    bool backFaceCulling(const Vector &, const Normal &);

// Calcul couleur pixel
    Vector calculeLumieresPhong( Point&, vector<Point> &, const PNTriangle& );

// Fonctions utilisées pour calculer la subdivision de phong
    Vector projecOrthoPhong(const Point &, const Point &, const Normal &);
    Normal phongNormalInterpolation(const double , const double , const PNTriangle &);
    Point phongTesselation(const double , const double , const PNTriangle &);

// Fonctions utilisées pour calculer la subdivision en PN triangles
    Normal calculNormaleUVPN(const Normal &n200, const Normal &n020, const Normal &n002,
                             const Normal &n110, const Normal &n011, const Normal &n101, const double u, const double v);
    double calculWij(const Point &P1, const Point &P2, const Vector &N);
    Point calculPointPN(const Point& P1, const Point& P2, const Vector &N);
    double calculVij(const Point &P1, const Point &P2, const Vector &N1,const Vector &N2);
    Normal calculNormalPN(const Point &P1, const Point &P2,  Vector &N1,const Vector &N2);

// Vérification point dans triangle
    bool memeCote(const Point &, const Point &, const Point &, const Point &);
    bool pointDansTriangle(const Triangle &, const Point &);

// Remplissage des triangles
    void sample( Image *, const PNTriangle& , const Transform& , const Transform& , vector<double> &);
    void remplissage1PointParPixel(Image *, const PNTriangle& , const Transform& , const Transform& , vector<double> &);
    void remplissage1PointParPixel2(Image *, const PNTriangle& , const Transform& , const Transform& , vector<double> &);

// Subdivisions
    void subdividePN   ( Image *, const PNTriangle& , const Transform&, const Transform&, vector<double> &, int);
    void subdividePhong( Image *, const PNTriangle& , const Transform&, const Transform&, vector<double> &, int);
    void subdivide     ( Image *, const PNTriangle& , const Transform&, const Transform&, vector<double> &);

    bool transform(Point&, const Point& , const Transform& , const Transform& );

/// ################################################################################### BACKFACE CULLING

bool backFaceCulling(const Vector &V, const Normal &N) { return ( Dot(Normalize(V), Normalize(N)) >= -0.7);}

/// ################################################################################### CALCUL DE COULEUR

Vector calculeLumieresPhong( Point&q, vector<Point> &l, const PNTriangle& triangle )
{
    Vector N = Vector(Normalize(triangle.getNormal()));
    Vector V(0, 0, 1);
    Vector Ia = intensiteAmbiante    * fMatiere;

    Vector intensiteTotale = Ia;
    for(unsigned int i = 0; i < l.size(); i++)
    {
        Vector L = Normalize(Vector(q-l[i]));
        double alpha = max(0.f, Dot(L,N));
        Vector R = 2 * alpha * N - L;

        double fMatiereSpeculaire = ((coeffSpeculaire +1)/M_PI)*cos(pow(acos(alpha),coeffSpeculaire));

        Vector Id = intensiteDiffuse     * fMatiere * alpha;
        Vector Is = intensiteSpeculaire  * fMatiereSpeculaire * pow(Dot(R,V),coeffSpeculaire);
        intensiteTotale += Id + Is;
    }
    for(int i = 0; i < 3; i++) if (intensiteTotale[i] > 255) intensiteTotale[i] = 255;

    return Vector(intensiteTotale);
}
/// ################################################################################### FONCTIONS UTILISEES POUR LA SUBDIVISIONS DE PHONG

Vector projecOrthoPhong(const Point &q, const Point &pi, const Normal &ni)
{
    Vector temp = Dot(Vector(q)-Vector(pi),Vector(ni)) * Vector(ni);
    return (Vector(q)-temp);
}

Normal phongNormalInterpolation(const double u, const double v, const PNTriangle &t)
{
    double w = 1 - u - v;
    return Normalize(Normal(Vector(u * t.getNormal().x, v * t.getNormal().y, w * t.getNormal().z)));
}

Point phongTesselation(const double u, const double v, const PNTriangle &t)
{
    double w = 1 - u -v;
    return (pow(u,2) * t.a + pow(v,2)* t.b + pow(w,2)* t.c +
            u*v*(projecOrthoPhong(t.b, t.a, t.na) + projecOrthoPhong(t.a, t.b, t.nb)) +
            v*w*(projecOrthoPhong(t.c, t.b, t.nb) + projecOrthoPhong(t.b, t.c, t.nc)) +
            w*u*(projecOrthoPhong(t.a, t.c, t.nc) + projecOrthoPhong(t.c, t.a, t.na)));
}




/// ################################################################################### FONCTIONS UTILISEES POUR LA SUBDIVISIONS EN PN TRIANGLES

Normal calculNormaleUVPN(const Normal &n200, const Normal &n020, const Normal &n002,
                         const Normal &n110, const Normal &n011, const Normal &n101, const double u, const double v)
{
    double w = 1 - u - v;
    return n200 * pow(w,2) + n020* pow(u,2) + n002 * pow(v,2) + n110 * w * u + n011 * u * v + n101 * w * v;
}

double calculWij(const Point &P1, const Point &P2, const Vector &N){return Dot(P2-P1,N);}

Point calculPointPN(const Point& P1, const Point& P2, const Vector &N)
{
    Point temp(N * calculWij(P1,P2,N) );
    return Point(((2 * P1 + P2) - temp)/3);
}

double calculVij(const Point &P1, const Point &P2, const Vector &N1,const Vector &N2) { return 2 * Dot(P2 - P1, N1 + N2)
                                                                                                  /Dot(P2 - P1, P2 - P1); }

Normal calculNormalPN(const Point &P1, const Point &P2,  Vector &N1,const Vector &N2)
{
    Vector temp(calculVij(P1,P2,N1,N2) * (P2 - P1));
    Point temp2(N1+N2 - temp);
    return Normal(Vector(temp2/ sqrt(pow(temp2.x,2)+pow(temp2.y,2)+pow(temp2.z,2)) ));
}
/// ################################################################################### CACULS SI UN POINT EST DANS UN TRIANGLE

bool memeCote(const Point &p1, const Point &p2, const Point &a, const Point &b)
{
  return (Dot(Cross(Vector(b-a),Vector(p1-a)), Cross(Vector(b-a),Vector(p2-a))) >= 0);
}

bool pointDansTriangle(const Triangle &t, const Point &p)
{
    return (memeCote(p, t.a, t.b, t.c) && memeCote(p, t.b, t.a, t.c) && memeCote(p, t.c, t.a ,t.b));
}

/// ################################################################################### APPLICATION DES TRANSFORMATIONS A UN POINT

bool transform( gk::Point& q, const gk::Point& p, const gk::Transform& mvp, const gk::Transform& viewport )
{
    gk::HPoint h;
    // projette le point
    mvp(p, h);

    // verifie que le point est visible
    if(h.isCulled())
        return false;

    // si le point est visible, finir la transformation
    q= viewport(h.project());
    return true;
}

/// ################################################################################### FONCTIONS DE REMPLISSAGE D'UN TRIANGLE

void sample( gk::Image *image, const gk::PNTriangle& triangle, const gk::Transform& mvp, const gk::Transform& viewport, std::vector<double> &zbuffer)
{
    // genere 20 points aleatoirement a l'interieur du triangle
    for(int i= 0; i < 100; i++)
    {
        // genere un point aleatoire a l'interieur du triangle
        gk::Point p;
        triangle.sampleUniform(sampler, p);

        // transforme le point
        gk::Point q;
        if(transform(q, p, mvp, viewport))
        {
             int cX = (int)q.y * largeur + (int)q.x;
             if( zbuffer[cX] > q.z)
             {
                zbuffer[cX] = q.z;

                image->setPixel( q.x, q.y, gk::Pixel(colorB, colorG, colorR) );
                nbPixelEcrit++;
             }
             else nbPixelNonEcrit ++;
        }
    }
    colorB = (int)(triangle.getNormal().x * 255)%255;
    colorG = (int)(triangle.getNormal().y * 255)%255;
    colorR = (int)(triangle.getNormal().z * 255)%255;
}



void remplissage1PointParPixel(gk::Image *image, const gk::PNTriangle& triangle, const gk::Transform& mvp , const gk::Transform& viewport, std::vector<double> &zbuffer)
{
    int nbPX, nbPY;
    gk::HPoint h1,h2,h3;
    Point p1,p2,p3, np;

    Triangle abc(triangle.a, triangle.b, triangle.c);

    BBox bbox = abc.getBBox();

    nbPX = sqrt(pow( bbox[0].x - bbox[1].x, 2)) /taillePixel;
    nbPY = sqrt(pow( bbox[0].y - bbox[1].y, 2)) / taillePixel;

    Vector AB(abc.b - abc.a);
    Vector AC(abc.c - abc.a);

    Vector pV = gk::Cross(AB,AC);

    if(nbPX >= 0 && nbPY >= 0)
    {
        for (int i = 0; i <= nbPX; i++)
        {
            for( int j = 0; j <= nbPY; j++)
            {
               double x = bbox[0].x + i *taillePixel;
               double y = bbox[0].y + j *taillePixel;
               double z =  ((-pV.x*(x - abc.a.x) - pV.y *(y - abc.a.y)) / pV.z) + abc.a.z;

                gk::Point p (x, y, z);

                if(pointDansTriangle(abc, p))
                {
                     gk::Point q;


                    if(transform(q, p, mvp, viewport))
                    {


                         int cX = (int)q.y * largeur + (int)q.x;

                         if( zbuffer[cX] > q.z) /// vérification Z- buffer
                         {
                            zbuffer[cX] = q.z;

                            /// calcul de la couleur du point à afficher
                            vector<Point> sourcesLumiereTransformees;
                            sourcesLumiereTransformees.reserve(sourcesLumiere.size());
                            for (unsigned int k = 0; k < sourcesLumiere.size(); k++ )
                            {
                                Point tempL;
                                transform(tempL, sourcesLumiere[k], mvp, viewport);
                                sourcesLumiereTransformees.push_back(tempL);
                            }

                            Vector couleur = calculeLumieresPhong(q, sourcesLumiereTransformees, triangle);

                            image->setPixel( q.x, q.y, gk::Pixel(couleur[0], couleur[1],couleur[2]) );
//                          image->setPixel( q.x, q.y, gk::Pixel(150,150,150 ));
                          image->setPixel( q.x, q.y, gk::Pixel(colorB, colorG, colorR) );
                            nbPixelEcrit++;
                         }
                    }
                }
            }
        }
    }

    colorB = (int)(abc.getNormal().x * 255)%255;
    colorG = (int)(abc.getNormal().y * 255)%255;
    colorR = (int)(abc.getNormal().z * 255)%255;

}

/// NON UTILISE (bug sur les jointures des triangles)
void remplissage1PointParPixel2(gk::Image *image, const gk::PNTriangle& triangle, const gk::Transform& mvp , const gk::Transform& viewport, std::vector<double> &zbuffer)
{
    Point aT, bT, cT;
    transform(aT,  triangle.a, mvp, viewport);
    transform(bT,  triangle.b, mvp, viewport);
    transform(cT,  triangle.c, mvp, viewport);
    Triangle abc =  Triangle(aT,bT,cT);

    BBox bbox = abc.getBBox();
    double dX = sqrt(pow( bbox[0].x - bbox[1].x, 2));
    double dY = sqrt(pow( bbox[0].y - bbox[1].y, 2));

    Vector AB(bT - aT);
    Vector AC(cT - aT);

    Vector pV = gk::Cross(AB,AC);

    vector<Point> grillePixel;

    double x = bbox[0].x, y = 0, z = 0;

    while( x < bbox[0].x + dX)
    {
        y = bbox[0].y;
        while(y < bbox[0].y + dY)
        {
            z =  ((-pV.x*(x - aT.x) - pV.y *(y - aT.y)) / pV.z) + aT.z;

                grillePixel.push_back(Point(x,y,z));
            y++;
        }
        x++;
    }


    for(unsigned int i = 0; i < grillePixel.size(); i++)
    {
        if(pointDansTriangle(abc, grillePixel[i]))
        {

             int cX = (int)grillePixel[i].y * largeur + (int)grillePixel[i].x;

             if( zbuffer[cX] > grillePixel[i].z) /// vérification Z- buffer
             {
                zbuffer[cX] = grillePixel[i].z;

                /// calcul de la couleur du point à afficher
                vector<Point> sourcesLumiereTransformees;
                sourcesLumiereTransformees.reserve(sourcesLumiere.size());
                for (unsigned int k = 0; k < sourcesLumiere.size(); k++ )
                {
                    Point tempL;
                    transform(tempL, sourcesLumiere[k], mvp, viewport);
                    sourcesLumiereTransformees.push_back(tempL);
                }

                Vector couleur = calculeLumieresPhong(grillePixel[i], sourcesLumiereTransformees, triangle);
                image->setPixel( grillePixel[i].x, grillePixel[i].y, gk::Pixel(couleur[0], couleur[1],couleur[2]) );
                nbPixelEcrit++;
             }
        }
    }
}


/// ################################################################################### FIN REMPLISSAGE



void subdividePN( gk::Image *image, const gk::PNTriangle& t, const gk::Transform& mvp , const gk::Transform& viewport, std::vector<double> &zbuffer, int nbPasses)
{
    Point pRa,pRb,pRc;

    // Définitions des points et normales des 3 points sommets du triangle
        Point P1(t.a); Vector N1(t.na); // b300 = P1
        Point P2(t.b); Vector N2(t.nb); // b030 = P2
        Point P3(t.c); Vector N3(t.nc); // b003 = P3

        Point b300(t.a);
        Point b030(t.b);
        Point b003(t.c);

    // Calculs des points de subdivision
        Point b210 = calculPointPN(P1,P2,N1);
        Point b120 = calculPointPN(P2,P1,N2);
        Point b021 = calculPointPN(P2,P3,N2);
        Point b012 = calculPointPN(P3,P2,N3);
        Point b102 = calculPointPN(P3,P1,N3);
        Point b201 = calculPointPN(P1,P3,N1);

        Point E = (b210 + b120 + b021 + b012 + b102 + b201) / 6;
        Point V = (P1+P2+P3) / 3;
        Point b111 = E +(E-V)/2;

    // Calculs des normales
        Normal n110 = calculNormalPN(P1,P2,N1,N2);
        Normal n011 = calculNormalPN(P2,P3,N2,N3);
        Normal n101 = calculNormalPN(P3,P1,N3,N1);
        Normal ncentre = calculNormaleUVPN(Normal(N1),Normal(N2),Normal(N3),Normal(n110),Normal(n011),Normal(n101), (1.f/3), (1.f/3));

    // Définitions des nouveaux PNTriangle de subdivision
        PNTriangle PN1(b300  , Normal(N1), b210, n110       , b201, n101);
        PNTriangle PN2(b210  , n110      , b120, n110       , b111, ncentre);
        PNTriangle PN3(b120  , n110      , b030, Normal(N2) , b021, n011);
        PNTriangle PN4(b120  , n110      , b021, n011       , b111, ncentre);
        PNTriangle PN5(b210  , n110      , b111, ncentre    , b201, n101);
        PNTriangle PN6(b201  , n101      , b111, ncentre    , b102, n101);
        PNTriangle PN7(b102  , n101      , b111, ncentre    , b012, n011);
        PNTriangle PN8(b012  , n011      , b111, ncentre    , b021, n011);
        PNTriangle PN9(b102  , n101      , b012, n011       , b003, Normal(N3));


    Vector V1(0, 0, 10); /// vecteur orienté vers l'observateur

    if(backFaceCulling(V1, t.getNormal()))
    {
         if( ((t.a.x != t.b.x) || (t.b.x !=t.c.x)) || ((t.a.y != t.b.y) || (t.b.y !=t.c.y)) || ((t.a.z != t.b.z) || (t.b.z != t.c.z)) )
            {

                if( transform(pRa,t.a,mvp,viewport) && transform(pRb,t.b,mvp,viewport) && transform(pRc,t.c,mvp,viewport))
                {
                    if(nbPasses == 0 || t.getArea() < 0.1)
                    {
                        remplissage1PointParPixel(image,PN1,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN2,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN3,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN4,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN5,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN6,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN7,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN8,mvp,viewport,zbuffer);
                        remplissage1PointParPixel(image,PN9,mvp,viewport,zbuffer);
                    }
                    else
                    {
                        nbPasses--;
                        subdividePN(image, PN1, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN2, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN3, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN4, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN5, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN6, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN7, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN8, mvp, viewport,zbuffer,nbPasses);
                        subdividePN(image, PN9, mvp, viewport,zbuffer,nbPasses);
                    }
                    nbTriangle += 9;
                }

                else if(t.getArea() > 0.01) /// on arrête si la surface devient trop petite
                {

                      /// on subdivise tant qu'un point du triangle est visible
                      if( transform(pRa,t.a,mvp,viewport) || transform(pRb,t.b,mvp,viewport) || transform(pRc,t.c,mvp,viewport))
                      {
                        subdividePN(image, PN1, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN2, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN3, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN4, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN5, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN6, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN7, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN8, mvp, viewport,zbuffer,0);
                        subdividePN(image, PN9, mvp, viewport,zbuffer,0);
                      }
                }
        }
    }
}


void subdividePhong( gk::Image *image, const gk::PNTriangle& t, const gk::Transform& mvp , const gk::Transform& viewport, std::vector<double> &zbuffer, int nombrePassesPhong)
{
    Point pRa,pRb,pRc;

    Point P1(t.a); Vector N1(t.na); // b300 = P1
    Point P2(t.b); Vector N2(t.nb); // b030 = P2
    Point P3(t.c); Vector N3(t.nc); // b003 = P3

    Point PM = phongTesselation(1.f/3,1.f/3,t);
    Normal nPM = phongNormalInterpolation(1.f/3,1.f/3,t);

    PNTriangle PN1(P1, Normal(N1), P2, Normal(N2), PM, nPM);
    PNTriangle PN2(P2, Normal(N2), P3, Normal(N3), PM, nPM);
    PNTriangle PN3(P3, Normal(N3), P1, Normal(N1), PM, nPM);

    Vector V1(0, 0, 10); /// vecteur orienté vers l'observateur

    if(backFaceCulling(V1, t.getNormal()))
    {
        if( ((t.a.x != t.b.x) || (t.b.x !=t.c.x)) || ((t.a.y != t.b.y) || (t.b.y !=t.c.y)) || ((t.a.z != t.b.z) || (t.b.z != t.c.z)) )
        {
            if( transform(pRa,t.a,mvp,viewport) && transform(pRb,t.b,mvp,viewport) && transform(pRc,t.c,mvp,viewport))
            {
                if(nombrePassesPhong == 0 || t.getArea() < 0.3)
                {
                    remplissage1PointParPixel(image,PN1,mvp,viewport,zbuffer);
                    remplissage1PointParPixel(image,PN2,mvp,viewport,zbuffer);
                    remplissage1PointParPixel(image,PN3,mvp,viewport,zbuffer);
                }
                else
                {

                    subdividePhong(image, PN1, mvp, viewport,zbuffer,nombrePassesPhong-1);
                    subdividePhong(image, PN2, mvp, viewport,zbuffer,nombrePassesPhong-1);
                    subdividePhong(image, PN3, mvp, viewport,zbuffer,nombrePassesPhong-1);
                }

                nbTriangle += 3;
            }

            else if(t.getArea() > 0.01) /// on arrête si la surface devient trop petite
            {

              /// on subdivise tant qu'un point du triangle est visible
              if( transform(pRa,t.a,mvp,viewport) || transform(pRb,t.b,mvp,viewport) || transform(pRc,t.c,mvp,viewport))
              {
                    subdividePhong(image, PN1, mvp, viewport,zbuffer,nombrePassesPhong-1);
                    subdividePhong(image, PN2, mvp, viewport,zbuffer,nombrePassesPhong-1);
                    subdividePhong(image, PN3, mvp, viewport,zbuffer,nombrePassesPhong-1);
              }
            }
        }
    }
}

void subdivide( gk::Image *image, const gk::PNTriangle& t, const gk::Transform& mvp , const gk::Transform& viewport, std::vector<double> &zbuffer)
{
    Point pRa,pRb,pRc;
    Vector V1(0, 0, 10); /// vecteur orienté vers l'observateur

    if(backFaceCulling(V1, t.getNormal()))
    {
         if( ((t.a.x != t.b.x) || (t.b.x !=t.c.x)) || ((t.a.y != t.b.y) || (t.b.y !=t.c.y)) || ((t.a.z != t.b.z) || (t.b.z != t.c.z)) )
            {
                nbTriangle++;
                if( transform(pRa,t.a,mvp,viewport) && transform(pRb,t.b,mvp,viewport) && transform(pRc,t.c,mvp,viewport))
                {
                     remplissage1PointParPixel(image,t,mvp,viewport,zbuffer);
                }
                else if(t.getArea() > 0.01) /// on arrête si la surface devient trop petite
                {
                     /// on subdivise tant qu'un point du triangle est visible
                      if( transform(pRa,t.a,mvp,viewport) || transform(pRb,t.b,mvp,viewport) || transform(pRc,t.c,mvp,viewport))
                      {
                            Point Mab = t.getUVPoint(0.5,  0);
                            Point Mbc = t.getUVPoint(0.5,  0.5);
                            Point Mca = t.getUVPoint(0,    0.5);
                            Normal NMab,NMbc,NMca;

                            NMab = (t.na + t.nb)/2; NMbc = (t.nb + t.nc)/2; NMca = (t.nc + t.na)/2;
                            PNTriangle nTriangle1(Mab,NMab,Mbc,NMbc,Mca,NMca);
                            PNTriangle nTriangle2(t.a,t.na,Mab,NMab,Mbc,NMbc);
                            PNTriangle nTriangle3(Mab,NMab,t.b,t.nb,t.c,t.nc);
                            PNTriangle nTriangle4(t.c,t.nc,Mab,NMab,t.a,t.na);

                            subdivide(image, nTriangle1, mvp, viewport,zbuffer);
                            subdivide(image, nTriangle2, mvp, viewport,zbuffer);
                            subdivide(image, nTriangle3, mvp, viewport,zbuffer);
                            subdivide(image, nTriangle4, mvp, viewport,zbuffer);
                      }
                }
            }
    }
}

int main( int argc, char **argv )
{

    // initialise le generateur de nombre aleatoire
    nbPixelEcrit = 0;
    nbPixelNonEcrit = 0;

    gk::Sampler::init(0);
    colorB = 255;
    colorG = 255;
    colorR = 0;
    nbDehors = 0;

    /// initialisation des sources de lumières
    sourcesLumiere.push_back(Point(15.f,0.f,100.f));
    sourcesLumiere.push_back(Point(10.f,2.f,100.f));

    std::vector<double> zbuffer;



    // cree l'image resultat
    gk::Image *image= new gk::Image(hauteur, largeur);

    // initialisation zbuffer
    zbuffer.assign(hauteur * largeur, 1);
    // charge un objet
    gk::Mesh *mesh= gk::MeshIO::read("bigguy.obj");
    if(mesh == NULL)
        return 1;


    // definit les transformations
    gk::Transform model;
    gk::Transform view       = gk::Translate( gk::Vector(0.f, 0.f, -50.f) );
    gk::Transform projection = gk::Perspective          (50.f, 1.f, 1.f, 1000.f);
    gk::Transform viewport   = gk::Viewport(image->width(), image->height());

    // compose la transformation complete
    gk::Transform mvp= projection * view * model;

    // recupere les triangles de la surface de l'objet
    int count= mesh->triangleCount();
    for(int i= 0; i < count; i++)
    {
        const gk::PNTriangle& triangle= mesh->getPNTriangle(i);

        /// choisir ici la subdivision
//        subdividePN(image, triangle, mvp, viewport, zbuffer,0);
        subdividePhong(image, triangle, mvp, viewport, zbuffer,1);
//        subdivide(image, triangle, mvp, viewport, zbuffer);
    }

    cout << "nbPixelEcrit   : "  << nbPixelEcrit << endl;
    cout << "triangle: "  << nbTriangle<< endl;

    // enregistre le resultat.
    gk::ImageIO::write(image, "output.bmp");
    return 0;
}
